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CALCULATION  OF  STATE  PROBABILITIES  FOR  A 


STOCHASTIC  LANCHESTER  COMBAT  MODEL 

by 

L.  Billard 

Florida  State  University 

1.  INTRODUCTION 

Lanchester  (1914)  presented  a  mathematical  model  describing 
the  conflict  situation  of  two  forces  in  combat  losing  units 
due  to  attrition.  This  model  was  deterministic  in  nature  and  as, 
for  example, in  Gye  and  Lewis  (1976)  this  approach  is  known  to 
break  down  in  certain  cases  most  especially  when  the  two  sides  are 
nearly  or  equally  matched.  In  these  circumstances,  it  is  more 
imperative  that  the  process  be  described  from  a  stochastic  out¬ 
look  . 

Therefore  Billard  (1979)  ..onsiders  how  stochastic  analogues 
of  many  deterministic  Lanchester-type  combat  models  can  be 
formulated.  Basically  the  processes  can  be  represented  as  either 
a  bivariate  death  process  or  a  bivariate  birth  and  death  process. 
That  paper  includes  a  brief  discussion  of  how  the  resultant 
Chapman-Kolmogorov  (differential-difference)  equations  can  be 
solved  to  provide  expressions  for  tl.e  underlying  state  prob¬ 
abilities.  The  techniques  employed  are  general  in  nature  and 
can  be  used  successfully  on  a  very  wide  ranging  array  of  model 


It  is  proposed  here  to  demonstrate  just  how  these 


techniques  can  be  employed.  This  is  done  using  the  original 
Lanchester  process,  a  process  of  great  interest  in  its  own 
right.  First,  however,  the  basic  mathematical  deterministic 
and  stochastic  models  for  the  Lanchester  process  is  reviewed 
briefly  in  Section  2.  A  discussion  of  the  Severo  (1969)  recur¬ 
sion  solution  technique  is  presented  in  Section  3  as  it  pertains 
to  the  special  nature  of  our  problem.  A  partitioning  scheme 
which  further  facilitates  exploitation  of  the  underlying 
structure  is  described  in  Section  4.  Finally,  in  Section  5 
the  solution  is  developed. 
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2.  THE  LANC HESTER  MODEL 


Let  B(t)  and  R(t)  be  the  number  of  units  (men,  tanks, 
etc.)  at  time  t  in  the  Blue  and  Red  force,  respectively.  Let 
particular  values  of  these  (integer  valued)  variables  be  b 
and  r,  respectively.  Let  B(0)  =  Bq  and  R(0)  =  Rq ,  assumed 
to  be  known  nonnegative  quantities.  The  two  forces  are  in  combat 
with  each  other  and  lose  units  by  attrition.  In  the  original 
Lanchester  (1914)  model  formulation,  the  force  sizes  at  time  t 
are  modelled  deterministically  according  to  the  equations 


db 

at  - 


-Ylr 


(1) 


and 


if  * 


f 


where  and  Y2  are  the  attrition  rates  of  the  Red  and  Blue 

forces,  respectively.  For  convenience, we  rescale  so  that  Yj  =  1 
and  y^  =  A  where  now  X  is  the  relative  effectiveness  of  the 
Blue  force  to  the  Red  force. 

A  stochastic  model  for  this  Lanchester  model  can  be  pre¬ 
sented  when  we  view  the  process  as  a  bivariate  pure  death 
process.  Thus  B(t)  and  R(t)  are  now  random  variables. 

Suppose 

pb,r(t)  =  Pr{B(t)  =  b'  =  r}- 
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Then,  the  corresponding  Chapman-Kolmogorov  (differential- 
difference)  equations  are 


4  Pb,r(t)  *  '  (Xr  +  b)  Pb,r(t)  +  XrPb+l,r(t>  +  bpb,r+l<t) 

for  (b,r)  €  A  *  { (b,r) :  0  £  b  £  Bq,  0  <  r  <  Rjj}.  For  boundary 
values  of  A,  certain  obvious  adjustments  are  necessary  in  (2a). 
when  b  =  BQ,  that  is,  for  points  (BQ,r)  in  A,  r  ^  0  and 
r  3*  Rq  the  Chapman-Kolmogorov  equations  are 


S  Vr<b> 


-(«  +  B0)Pgo<r(t)  +  B0pBoir+1(t) 


when  r  =  RQ,  that  is,  for  points  (b, Rq)  in  A,  b  f  BQ, 


*Pr,R0<fc)  * 


-(XRn  +  b)  pbR  (t)  +  XRnpKAl  D  (t)  ? 


O^b+1 , R, 


when  r  =  0,  that  is,  for  points  (b,0)  in  A, 


3t  Pb,0(t)  =  bPb,l(t) 


when  b  =  0,  that  is,  for  points  (0,b)  in  A, 


3t  PQ,r(t)  =  XrPl,r(t)  ; 


and  when  (b,r)  =  (Bq,Rq), 


5t  PB0,R0(t'  *  -(XR0  +  V  W0  ' 


*  \  * 


,  (2a) 


Thus , 


(2b) 


(2c) 


(2d) 


(2e) 


(2f ) 


These  equations  (2)  can  be  solved  by  adopting  the  tech¬ 
nique  of  Severo  (1969) .  The  first  step  then  is  to  identify 
each  point  (b,r)  in  A  by  a  counting  coordinate  k  given  by 

k  *  k (b,r)  =  (BQ  +  1)  (RQ  +  l )  -  fc(R0  +  1)  -  r  . 

Then  if  we  set 

yk(t)  *  Pb,r(t)  '  (3) 

we  may  rewrite  (2)  in  terms  of  the  new  coordinate  k.  We  suffice 
by  writing  (2a)  only,  the  other  parts  of  (2)  following  similarly. 
Thus,  equation  (2a)  becomes 

yk(t)  =  -Ur  +  b)  yk(t)  +  Ar  yk_R  _i  (t)  +  byk-1(t)  .  (4) 

Or,  in  matrix  terms 

^  £(t)  =  B^(t)  .  (5) 

Clearly,  from  (4),  the  matrix  of  coefficients  B  is  lower  tri¬ 
angular.  Therefore  we  can  use  Severo' s  recursive  result  to  give 
us 

£(t)=  Ce(t)  (6) 

where  e(t)  is  the  (BQ  +  1) (RQ  +  1)  *  1  vector  with  elements 
exptb^t)  with  b^  being  the  kth  diagonal  element  of  B.  The 
elements  c(i,j),  i,  j  =  !»...»  (Bq  +  1) (Rq  +  1) »  of  the  matrix 
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C  can  be  found  from  Severo  (1969,  Theorem  1) .  That  is, 


c (i, j )  = 


al  ' 

i-1 

a.  -  T  c (i,u) 
1  u*l 


i  -  j  -  1. 


i  =  j  >  2, 


b’ (i, i-1)  C(i-l,j)  h ( j , i-1)  ,  i  >  j, 


0  , 


i  < 


where 


c(i,j)  =  cQ(i,  j)  +  +•••+  (i,  j)  ; 

i  >  j  , 


(i' j)  "  (bil'  • • *  '  bij)  ' 


(7) 


(8) 


with  b^  the  (ij)th  element  of  B  ; 


h’(j,i)  =  (6  0  (b  j  -  b^,...,  5  i_j_i  *  i  >  j  , 


with 


(9) 


&L(y)  = 


(i!/yi+1)  l  (-1)1"3  (yt)Vj!  ,  y  /  0  , 


j=0 


ti+1/(i+l) , 


y  =  0; 


(C (i, j) ) 


rs 


cs-l(r'j)  ' 
0 


r  1  j  +  s  -  1 , 
otherwise  , 


for  r  *  1,...,  i  and  s  =  1,...,  i-j+1  ;  and  ak  =  yk<0), 
k  =  1, . . .,  (Bq+1) (R0+l) • 


(10) 


*  **».*•  v  <v  j  v  * 


3.  THE  SOLUTION  MATRIX  C 


The  expressions  for  finding  C  given  in  the  previous 
section  are  the  general  ones  given  by  Severo.  However,  a  close 
examination  of  our  particular  problem  shows  that  some  interest¬ 
ing  simplifications  can  be  obtained. 

Typically,  a  =  (1,0, This  corresponds  to 
B(0)  =  Bq  and  R(0)  =  Rq  with  probability  one.  We  shall  assume 
this  holds  in  the  rest  of  this  paper.  The  adjustment  for  the 
general  a  case  is  trivial. 

We  notice  that  b(i,i-l)  is  simply  the  ith  row  of  the 
subdiagonal  matrix  of  B.  From  (2)  (or  (5)),  in  our  case  this 
vector  has  at  most  two  "non-zero"  elements  corresponding  to  the 
coefficients  of  y^^tr)  an<^  ^k-R  -l^'  Note  that  in  some 
instances  these  values  may  "happen  to  be"  zero.  All  other 
elements  of  b(i,i-l)  are  zero.  We  exploit  this  fact  in  the 
actual  determination  of  the  elements  of  C. 

A  further  simplification  occurs  when  the  coefficients 
of  B  are  time  independent  as  is  the  situation  for  our  Lanchester 
model.  We  first  note  that  the  recursiveness  in  the  solution 
enters  through  the  C  (i-l,j)  matrix.  We  also  observe  that 
in  general  we  can  express  obtained  elements  c(i,j)  as  a  poly¬ 
nomial  in  t,  viz., 

c(i,j)  =  cQ(i,j)  +  tc  ^  ( i ,  j )  +  •••  +  t1  ^c^jtifj)  . 
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present  model,  the  elements  c(i,j)  are  time  independent.  That 
is,  all  elements  in  columns  other  than  the  first  are  zero.  Hence, 
in  h(i-l,j)  we  need  only  determine  the  first  element  6g(bj-b^). 
As  a  final  comment,  if  b^  =  b^,  then  6g(0)  is  a  function  of 
t.  However,  typically  this  situation  arises  when  we  are  in  a 
row  i  whose  b(i,i-l)  C  (i-l,j)  elements  are  all  zero.  Thus, 
the  independence  of  c(i,j)  on  t  remains  intact. 

By  exploiting  these  features,  we  can  determine  more 
readily  the  elements  of  C.  Before  doing  this  (in  Section  5), 
we  shall  present  a  partitioned  form  of  B  and  hence  C 
which  allows  us  to  observe  much  of  this  underlying  structure  at 
a  glance. 
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4.  A  PARTITIONING  SCHEME 

First  we  note  that  the  counting  coordinate  k  orders 
the  points  (b,r)  in  A  as  follows:  (B0,RQ),  (Bq,R0-1) , . . . , (BQ,0) ; 
(Bq-1 , R q ) # • • » i  (Bq— 1 / 0) ; . . . ;  (0/ Rq ) ,  •  • . ,  (0/0) .  Thus  we  partition 
these  into  groups  whose  b  coordinate  is  common.  That  is,  we 
may  write 

B  =  (b^)  ,  u,  v  =  1,...,  Bq+1, 

and  the  "elements,"  or  submatrices,  byv  have  elements 

buV  =  (buv(p,q)),  p,  q  =  1,...,  RQ+1  . 

This  is  illustrated  in  Table  1  for  the  case  (Bq,Rq)  =  (4,2). 
Thus,  the  k  =  k(b,r)th  row  of  B  corresponds  to  the 
p  =  (Rg-r+Dth  row  of  the  u  =  (Bg-b+Dtli  submatrix  of  B; 
similarly  for  the  columns  of  B. 

To  insert  the  appropriate  coefficients  in  B,  we  refer 
to  equation  (2)  (or  (4)).  It  is  immediately  apparent  that  the 
coefficients  of  p^  r(t)  (or  y^ ( t) )  constitute  the  diagonal 
elements  of  B;  the  coefficients  of  p^  r+^(t)  (or  y^^tt)) 
consistute  the  off-diagonal  elements  of  the  diagonal  submatrices; 
the  coefficients  of  pfa+1  r(t)  (or  yk_R  _]_(t))  constitute 
the  diagonal  elements  of  the  off-diagonal  submatrices;  and  all 
other  elements  are  zero.  If  we  assert  that  of f-diagonal  elements 
of  a  submatrix  cannot  go  "outside"  that  submatrix,  then  the 
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adjustments  to  equation  (2a)  which  produced  the  remaining  equations 
of  (2)  when  dealing  with  boundary  values  of  (b,r)  are  auto¬ 
matically  taken  care  of.  Thus,  for  example,  when  (Bq,Rq)  ■  (4,2) 
the  equation  relating  to  the  boundary  point  (2,2)  does  not  have 
a  term  corresponding  to  p2  3(t)  (or  y^^tt)).  so*  our  ® 

matrix  would  have  an  entry  in  the  spot  indicated  by  *  in 
Table  1 . 

Thus,  in  our  Lanchester  model,  B  has  elements 


buu(p'p) 


-U(R0-p+l)  +  (Bq-U+1)  } 


u  1,..., Bq  ,  p = 1 , • . • , Rq 


0, 


u  =  Bq+1  and/or  p  =  Rq+1 


buu (p» P“ 1 )  =  (®q”^^1^ »  u  -  1, . . . ,Bq+1,  p—  2,..., Rq+1 { 

and 

bu  ^^^(p,p)  A (Rq~p+1 ) ,  u  2,..., Bq+1,  p  1,..., Rq  +1 . 

Clearly  then  such  a  partitioning  scheme  facilitates  the 
construction  of  the  matrix  of  coefficients  B.  We  shall  see 
that  the  analogous  partitioning  of  the  solution  matrix  C 
also  facilitates  the  derivation  of  its  elements.  Thus,  we  write 

-  =  (cuv)  '  u,  v  »  1,...,Bq+1, 

and  the  "elements,"  or  submatrices,  c 


have  elements 


5.  CALCULATION  OF  C 


Quite  clearly,  from  (6),  once  the  matrix  C  has  been 
determined,  the  state  probabilities  follow  immediately.  As 
indicated  earlier  we  exploit  the  underlying  structure  of  B 
and  the  accompanying  partitioning  scheme.  For  illustrative 
purposes,  let  us  take  the  particular  case  (BQ,R0)  *  (4,2) 
and  A  =  .8.  The  corresponding  B  matrix  of  coefficients  is 
given  in  Table  1 . 

From  equation  (7) ,  it  is  clear  that 

cuv(p,q)  =  0 

whenever  u  <  v  and/or  p  <  p.  That  is,  the  upper  diagonal 
submatrices  cuv  are  zero  and  the  upper  diagonals  of  the 
remaining  submatrices  are  also  zero. 

Let  us  take  the  elements  in  the  first  column  for 
(v,q)  =  (1,1) . 

First,  from  the  initial  condition,  we  have 

Cj^(l,l)  =  =  1  • 

Let  us  proceed  with  the  (1,1)  elements  of  the  submatrices 
(u,l),  u  =  1,...,  Bq+1*  We  find  from  Severo, 


(12) 


where  "-"  indicates  the  presence  of  a  "nonzero"  element  (which 
in  a  particular  case  may  happen  to  equal  zero) .  However  because 
of  the  zeros  in  b(i,i-l),  or  in  c  (i-l,j),  we  do  not  need  to 
determine  these  quantities.  Thus,  equation  (12)  reduces  to 

c21(l,l)  =  (1.6)  c11(l,l)/(b11(l,l)  -b22(l,l)} 

*  -1.6  . 


In  general, 

cul(l,l)  =  (1.6)  cu_1Jc(l,l)/{b11(l,l)  -buu(l,l)}  ; 

or, 

cul(l,l)  =  (-1) U_1 (1 . 6 ) u-1/u !  ,  u  =  1,...,  Bq+1  . 

Let  us  consider  the  column  (v,q)  =  (2,1).  Now, 
3 

c22(l,l)  =  -  l  c21(l,q)  -  1.6  +  0  ♦  0  -  1.6  . 


13 


Proceeding  as  before,  we  have 


=  (1.6)  c22(l,l)/{b22(l,l)  -  b33(l,l)} 


=  -2.56  . 


In  general, 

cu2(l,l)  =  (1.6)  cu_1  2(l,l)/{b22(l,l)  -  buu(l,l)} 

=  (-l)u  (1.6) U-1/ (u-1) !  ,  u=  2 , . . . , Bq+1 . 

Continuing  in  this  manner,  we  can  calculate  all  the 
cuv(l,l)  elements.  The  essential  feature  to  note  here  is 
that  in  order  to  determine  a  particular  cuv(l,l)  element, 
first  glance  at  the  Severo  result  equation  (7)  suggests  that 
all  previous  row  elements  in  that  column  of  C  are 
necessary,  whereas  by  exploiting  the  special  nature  of  our 
problem  we  only  require  the  previous  cu_j.  v(l»l)  element. 

Let  us  now  consider  the  elements  cuv(2,l).  In 
particular,  let  us  take  the  first  column  (v,q)  =  (1,1). 
Then,  from  (7), 
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c11(2,l)  -  4cnU,l)/{bu<l,l)  -  blx(2,2)} 
=  -5.0. 


Next, 


/  0 

-5.0 

c21(2,1)  =  (0  .80  3) 


^  /  5rt(*)\ 


\” 


1.6  0  0  0 


l  •/ 


{(.8)  c11(2,l)  +  3c, ,  (lajj/tb,,  (1,1)  -  b„<2,2)} 


'21 


11' 


'22 


=  4.889. 


In  general. 


cul * 2 ' 1 ^ 


=  {(.8)  cu-11(2,l)  +  (Bn-u+l)  ctll  (l,l)}/{b11  (1,1)  -  b„„(2,2))  , 


ul 


'11 


uu 


u  —  1,...,  Bq+1 


Similarly,  we  can  determine  that 


Cu2{2'1) 


{(‘8)  cu-1,2(2,1)  +  <Bo~u+1)  cu2(l,l)}/{b22(l,l)  -buu(2,2)}  , 


u  —  2,...,  Bq+1  . 
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Likewise  we  can  calculate  all  the  cuv(2,l)  elements.  Again 
an  essential  feature  is  that  to  find  a  particular  cuv(2,l) 
we  need  only  the  two  previous  elements  cu-1  v ( 2  # 1)  and 
cUu<l,l>  and  not  all  the  previous  elements  in  that  column. 

Therefore  by  repeating  this  process,  we  can  calculate 
the  matrix  C.  That  is,  we  determine  the  diagonal  element 

cuu(p,p)  -  -  cuv(p,q)  -  ^  =uu(p.q) 

Once  each  diagonal  element  is  calculated,  the  other  elements  in 
that  same  column  (v,q)  »  (u,p)  can  be  determined  according  to 

cuv(p'q) 

=  {X(RQ-p+l)  cu_1>v(p,q)  +  (B0-u+l)  c^tp-l^D/Cb^lqrql-b^lp.p)) 

(13) 

for  u  =  v,...,  Bq+1,  v  =  1,...,  Bq+1 ,  p  =  q, . . . ,R»+1, 
q  =  1,...,Rq+1,  and  where  equation  (11)  applies  when  applicable 
in  (13)  . 

For  the  computer  user  what  would  be  a  potential  storage 
problem  if  all  the  previous  cuv(p,q)  elements  were  needed, 
does  not  cause  any  undue  concern  since  only  the  two  appropriate 
values  need  be  retained.  This  assumes  the  program  has  been 
set  up  to  print  out  the  values  of  C  as  they  are  determined 
and  that  a  running  sum  of  the  row  elements  is  retained  (for  the 
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calculation  of  subsequent  diagonal  elements;  see  equation  (7) ) 
The  complete  matrix  C  for  our  example  is  shown  in 

Table  2. 
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6.  THE  STATE  PROBABILITIES 

Once  the  elements  of  C  have  been  calculated,  the 
state  probabilities  r<t)  follow  readily  from  (3)  and  (6). 
Thus,  for  example, 

Pr  (B(t)  *  2,  R(t)  =  1}=  P2,i<t> 

*  -2.311e"5,6t+  1.6e“4,8t+  5.511e"4*6t-  4.089e’3*8t 
-  3.200e”3*66  +  2 .489e”2  *8t  , 

and 

Pr  {Red  forces  lose  no  units} 

4 

=  l  Pr  {B  (t)  =  b,  R(t)  =  2} 
b*0 

=  0.192e”5'6t  +  0.376e‘4,6t  +  0.142e"3*6t  +  0.263e"2,6t  +  0.027 
In  general, 

PrtB(t)  »  b,  R(t)  »  r) 

Vb  B0+1 

=  l  CB0-b+l,V(Vr+1'q)  exP(tbvv(q,q), 

R0-r+1 

+  *  CBn-b+l,Bft-b+l(Vr+1'q)  eXP(tbBn-b+l,Bn-b+l(q»q>)‘ 

q«l  0  0  0  0 

Or,  more  simply,  in  terms  of  the  k  notation  used  before 
partitioning. 
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Pb,r<t>  5  *k(t) 


k 

l  c(k,j)  exp(b.t)  . 
j«l  3 


A  knowledge  of  the  state  probabilities  permits  the 
derivation  of  other  quantities  of  interest,  for  example,  the 
expected  force  sizes. 
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